#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "sign.h"
#include "materialfunction.h"

double alpha,m,n;
double theta_s = 1.0;
double theta_r = 0.0025;
double moisturecontent(double material[][3], double H,int mcase){
    double moi = 0.0;

    switch (mcase){
    case 1:
        alpha   = material[0][0];
        m       = material[0][1];
        n       = material[0][2];
        moi     = (theta_s - theta_r)/pow((pow(fabs(alpha*H),n) + 1),m )+theta_r;
        break;
    case 2:
        alpha   = material[1][0];
        m       = material[1][1];
        n       = material[1][2];
        moi     = (theta_s - theta_r)/pow((pow(fabs(alpha*H),n) + 1),m )+theta_r;
        break;
    case 3:
        alpha   = material[2][0];
        m       = material[2][1];
        n       = material[2][2];
        moi     = (theta_s - theta_r)/pow((pow(fabs(alpha*H),n) + 1),m )+theta_r;
        break;
    default:
        break;
      }

        return moi;

}

double conductivity(double material[][3],double H,int mcase){
    double con = 0.0;
    switch (mcase){
    case 1:
        alpha   = material[0][0];
        m       = material[0][1];
        n       = material[0][2];

        con     = alpha*m*n*pow(fabs(alpha*H),n - 1)*signf(alpha*H)*(theta_r - theta_s)/pow((pow(fabs(alpha*H),n) + 1),m + 1);
        break;
    case 2:
        alpha   = material[1][0];
        m       = material[1][1];
        n       = material[1][2];
        con     = alpha*m*n*pow(fabs(alpha*H),n - 1)*signf(alpha*H)*(theta_r - theta_s)/pow((pow(fabs(alpha*H),n) + 1),m + 1);
        break;
    case 3:
        alpha   = material[2][0];
        m       = material[2][1];
        n       = material[2][2];
        con     = alpha*m*n*pow(fabs(alpha*H),n - 1)*signf(alpha*H)*(theta_r - theta_s)/pow((pow(fabs(alpha*H),n) + 1),m + 1);
        break;
    default:
        break;
    }
    return con;
}

double permeability(double material[][3],double k0[],  double H,int mcase){
    double per = 0.0;

    switch (mcase){
    case 1:
        per = k0[0]*pow(moisturecontent(material,H,1),4.0);
        break;
    case 2:
        per = k0[1]*pow(moisturecontent(material,H,2),4.0);
        break;
    case 3:
        //printf("k02 %f\n",k0[2]);
        //printf("phi02 %f\n",phi[2]);
        //printf("moi %e\n",pow(moisturecontent(material,H,3)/phi[2],4.0));
        per = k0[2]*pow(moisturecontent(material,H,3),4.0);
        //printf("k02 %e\n",per);
        break;
    default:
        break;
    }
    return per;
}
